% Calculate the maximum values presented in Table 1.

clear
clc

% Designate working directories
dir = '.../Replication package/';
cd(dir);
addpath(genpath(strcat(dir,"Code/src")))
addpath(genpath(strcat(dir,"Datasets")))

% Set seed
seed = 2023 ;
rng(seed,'twister');

lambdaset = 0.3:-0.1:0;
a_grid = [1 2 3 4 5];
swarmsize = 5000;
theta = -0.028;
eps_constant = 1e-2;  % cutoff value for determining numerical zeros
pso_nmax = 100;       % maximum number of PSO trials for non-decreasing solutions w.r.t. λ

data = table2array(readtable('USAY.txt'));
data = rmmissing(data(:, 1:8));
T = size(data, 1) - 2;
Y = data(3:end, 6);
X = data(3:end, 8);
Z = [data(1:(end-2), 3:6)];
datau = Y - mean(Y) - theta*(X - mean(X));

max_vals = zeros(length(a_grid), length(lambdaset));

for i = 1:length(a_grid)
    for j = 1:length(lambdaset) 
        [stat,gamma,~,rtime] = pen_Bierens_test(datau,Z,a_grid(i),lambdaset(j),eps_constant,pso_nmax,'swarmsize',swarmsize,'dm',1);
        max_vals(i, j) = stat;
    end
end

disp(max_vals)

save(strcat(dir,"max_vals.mat"),"max_vals","lambdaset","a_grid","swarmsize")